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Abstract 

We consider the problem of globally minimizing the sum of many rational func- 
■ tions over a given compact semialgebraic set. The number of terms can be large 

(10 to 100), the degree of each term should be small (up to 10), and the number 
^ I of variables can be large (10 to 100) provided some kind of sparsity is present. We 

I^H ■ describe a formulation of the rational optimization problem as a generalized mo- 

Tlj" , ment problem and its hierarchy of convex semidefinite relaxations. Under some 

I conditions we prove that the sequence of optimal values converges to the globally 

optimal value. We show how public-domain software can be used to model and 
. solve such problems. 

(-H ■ Keywords: rational optimization; global optimization; semidefinite relaxations; sparsity. 

! AMS MSG 2010: 46N10, 65K05, 90C22, 90C26. 

^: 

1 Introduction 

' Consider the optimization problem 



TV 



r := inf J]/.(x) (1) 



O ■ i=l 

over the basic semi-algebraic set 

K := {x G M" : g,{^) > 0, j = 1, . . . , m }, (2) 

X ■ 

■ for given polynomials gj G IR[x], j = 1, . . . ,m, and where each term /j : 

- ■ ■' rational function 

X ^ /, X := — -, 



^Universite de Toulouse; Mines Albi; Institut Clement Ader (ICA); Campus Jarlard, F-81013 Albi, 
France 

^CNRS; LAAS; 7 avenue du colonel Roche, F-31077 Toulouse, France; Universite de Toulouse; UPS, 
INSA, INP, ISAE; LAAS; F-31077 Toulouse, France 

•^Faculty of Electrical Engineering, Czech Technical University in Prague, Technicka 4, CZ-16607 
Prague, Czech Republic 

^Institut de Mathematiques de Toulouse (IMT), Universite de Toulouse, France 

*The first author acknowledges support by the French National Research Agency (ANR) through 
COSINUS program (project ID4CS ANR-09-COSI-005). The second author acknowledges support by 
Research Program MSM6840770038 of the Czech Ministry of Education and Project 103/10/0628 of the 
Grant Agency of the Czech Republic. 



1 



with Pi, qi G M[x] and > on K, for each i = 1, . . . , A^. 

Problem ([1]) is a fractional programming problem of a rather general form. Nevertheless, 
we assume that the degree of each fi and gj is relatively small (up to 10), but the number 
of terms can be quite large (10 to 100). For dense data the number of variables n 
should also be small (up to 10). However, this number can be also quite large (10 to 100) 
provided that the problem data feature some kind of sparsity (to be specified later). Even 
though problem ([1]) is of self-interest, our initial motivation came from some applications 
in computer vision, where such problems are typical. These applications will be described 
elsewhere. 

In such a situation, fractional programming problem ([1]) is quite challenging. Indeed, 
we make no assumption on the polynomials pi, qi whereas even with a relatively small 
number of fractions and under convexity (resp. concavity) assumptions on pi (resp. g^), 
problem ([T]) is hard to solve (especially if one wants to compute the global minimum); see 
for example the survey pl)j and references therein. 

We are interested in solving problem ([T]) globally, in the sense that we do not content 
ourselves with a local optimum satisfying first order optimality conditions, as typically 
obtained with standard local optimization algorithms such as Newton's method or its 
variants. If problem ([1]) is too difficult to solve globally (because of ill-conditioning and/or 
too large a number of variables or terms in the objective function), we would like to have 
at least a valid lower bound on the global minimum, since upper bounds can be obtained 
with local optimization algorithms. 

One possible approach is to reduce all fractions Pi/qi to same denominator and obtain a 
single rational fraction to minimize. Then one may try to apply the hierarchy of semidef- 
inite programming (SDP) relaxations defined in [6], see also [9l Section 5.8]. But such 
a strategy is not appropriate because the degree of the common denominator is poten- 
tially large and even if n is small, one may not even implement the first relaxation of 
the hierarchy, due to the present limitations of SDP solvers. Moreover, in general this 
strategy also destroys potential sparsity patterns present in the original formulation ([T]), 
and so precludes from using an appropriate version (for the rational fraction case) of the 
sparse semidefinite relaxations introduced in [11] whose convergence was proved in [7] 
under some conditions on the sparsity pattern, see also P Sections 4.6 and 5.3.4]. 

Another possibility is to introduce additional variables r, (that we may call liftings) with 
associated constraints 

and solve the equivalent problem: 

N 

r := inf^ ^r, (3) 

which is now a polynomial optimization problem in the new variables (x, r) G IR^xM^, and 
where the new feasible set K = K x {(x, r) G M"^^ : rjgj(x) —pi{x.) > 0} is modeling the 
epigraphs of the rational terms. The sparsity pattern is preserved and if K is compact one 
may in general obtain upper and lower bounds fi,r^ on the rj so as to make K compact 
by adding the quadratic (redundant) constraints (r, — rj)(rj — r,) > 0, i = 1,...,N, 
and apply the sparse semidefinite relaxations. However, in doing so one introduces A^ 
additional variables, and this may have an impact on the overall performance, especially 
if A^ is large. In the sequel this approach is referred to as the epigraph approach. 



2 



The goal of the present paper is to circumvent all above difficulties in the following two 
situations: either n is relatively small, or n is potentially large but some sparsity is present, 
i.e., each /j and each gj in ([T]) is concerned with only a small subset of variables. In the 
approach that we propose, we do not need the epigraph liftings. The idea is to formulate 
([1]) as an equivalent infinite-dimensional linear problem which a particular instance of 
the generalized moment problem (GMP) as defined in [8], with N unknown measures 
(where each measure is associated with a fraction Pi/qi). In turn this problem can be 
easily modeled and solved with our public-domain software GloptiPoly 3 [5], a significant 
update of GloptiPoly 2 [1] . In the sequel this approach is referred to as the GMP approach. 

The outline of the paper is as follows. In Section H] we introduce the SDP relaxations ffist 
in the case that n is small and the data are dense polynomials. Then in Section Owe extend 
the SDP relaxations to the case that n is large but sparsity is present. In SectionlHwe show 
how the GMP formulation can be exploited to model the SDP relaxations of problem ([T]) 
easily with GloptiPoly 3. We also provide a collection of numerical experiments showing 
the relevance of our GMP approach, especially in comparison with the epigraph approach. 



2 Dense SDP relaxations 

In this section we assume that n, the number of variables in problem ([1]), is small, say up 
to 10. 



2.1 GMP formulation 

Consider the infinite dimensional linear problem 

N 

f 



inf / pidj^i 

s.t. / qidfjii = 1 (^4) 



x"g,d/ii = / x"girf/ii, Va G N", i = 2, . . . , AT, 

K JK 

where Ai(K) is the space of finite Borel measures supported on K. 

Theorem 2.1 Let K C M" m ^ be compact, and assume that qi > onJi, i = 1, . . . , N . 

Then f = f*. 

Proof: We ffist prove that f* > f. As / = J2iPi/li continuous on K, there exists a 
global minimizer x* G K with /(x*) = /*. Define fii := gj(x*)~^(5x*, i = 1, . . . ,N, where 
5x* is the Dirac measure at x*. Then obviously, the measures (/Xj), i = 1,...,N, are 
feasible for (jl]) with associated value 

2 = 1 1 = 1 

Conversely, let (fii) be a feasible solution of (jl]). For every i = 1, . . . , A^, let dui be the 
measure qidfii, i.e. 

iyi{B) := / gi(x)d/ii(x) 
JKnB 
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for all sets B in the Borel cr-algebra of M", and so the support of z/j is K. As measures on 
compact sets are moment determinate, the moments constraints of (jl]) imply that z/j = z/i, 
for every i = 2, . . . ,N, and from J^^ Qidfii = 1 we also deduce that ui is a probability 
measure on K. But then 

iv „ ^ r ^ r 

where we have used that / > /* on K and vi is a probability measure on K. □ 

We next make the following assumption meaning that set K admits an algebraic certificate 
of compactness. 

Assumption 2.1 The set K C M" m ^ is compact and the quadratic polynomial x i— )■ 
M — II can he written as 

m 

M - ||x|p = (To + ^crjgj, 
i=i 

for some polynomials aj G M[x], all sums of squares of polynomials. 



2.2 A hierarchy of dense SDP relaxations 

Let Yi = (yia) be a real sequence indexed in the canonical basis (x") of M[x], i = 1, . . . , N, 
and for every k eN, let N'^ := {a G : J2j < 

Define the moment matrix Mfc(yj) of order k, associated with y, whose entries indexed 
by multi-indices /3 (rows) and 7 (columns) read 

[Mk{y^)]p,^:=y.iP+^), V/3, 7eN^, 

and so are linear in yj. Similarly, given a polynomial ^'(x) = S'qx", define the localising 
matrix Mk{gyi) of order fc, associated with y and g, whose entries read 

[Mkigyi)]^,^ ■■= ^gayi{a+p+^), V/3,7 g N^. 

a 

In particular, matrix Mo{gyi) is identical to Ly-{g) where for every i, Ly. : R[x] — t- M is 
the linear functional defined by: 

9 ^ LyXg) := ^aVia, ^9 e M[x]. 

Let Ui := [(deggi)/2'l, i = 1,...,A^, rj := [(deg5'j)/2] , j = l,...,m, and with no loss 
of generality assume that Ui < U2 < ■ ■ ■ < u^. Consider the hierarchy of semidefinite 
programming (SDP) relaxations: 

N 

/,* = inf J2^yM 

s.t. Mk{yi)hO, t = l,...,N (5) 

Mk-r, i9jyi) ^0, i = l,...,N, j = l,...,m 

Ly,(x-g,) = Ly,(x-gi), Va G N^(,_„,^), ^ = 2, . . . , AT. 
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Theorem 2.2 Let Assumption \2.1\ hold and consider the hierarchy of SDP relaxations 
Then it follows that 

(a) fkt f* ask^oo. 

(b) Moreover, if (yf) is an optimal solution of and if 

rankMfc(yf ) = rankMfc_„^(yf) =: R, z = 1, . . . , iV 
then fl = f* and one may extract R global minimizers. 

Proof: The proof of (a) is classical. One first prove that if (yf) is a nearly optimal 
solution of ([5]), i.e. 

^ 1 

then there exists a subsequence [kp) and a sequence yj, i = 1, . . . , A^, such that 

\im yt = y^a, VaGN",z = l,...,iV. 

From this pointwise convergence it easily follows that for every i = 1,. . . ,N and j = 
1, . . . ,m, 

Mk{yi)hO, Mk{gjyi)hO, = 0, 1, . . . 

By Putinar's theorem ^ Theorem 2.14] this implies that the sequence y, has a repre- 
senting measure supported on K, i.e., there exists a finite Borel measure on K such 
that 

LyAf) = I fdfi^, V/GM[x]. 
Jk 

Moreover, still by pointwise convergence, 

Ly,(g,x") = / x'^g,(x)rf/i, = Ly,(gix") = / x'^gi(x)rf/ii, Va G N". (6) 
Jk Jk 

Therefore, let dui := qi{x.)dfii which is a probability measure supported on K. As K is 
compact, by dH]), z^j = z^i for every i = 1, . . . , N. Finally, again by pointwise convergence: 

f* > jim fl = lim Vl = V / pic//Xi 

i=l i=l "'^ 

= V / —q^df^i = —di^i 



which proves (a) because fl is monotone non- decreasing. In addition, vi is an optimal 
solution of (jl]) with optimal value /* = /. 

Statement (b) follows from the flat extension theorem of Curto and Fialkow [21 Theorem 
3.7] and each y^ has an atomic representing measure supported on R points of K. □ 
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3 Sparse SDP relaxations 



In this section we assume that n, the number of variables in problem ([T]), is large, say 
from 10 to 100, and moreover that some sparsity pattern is present in the polynomial 
data. 



3.1 GMP formulation 

Let Jo := {I, . . . ,n} = U^^^/j with possible overlaps, and let R[xfc : A; G Jj] denote the ring 
of polynomials in the variables x^, k & li. Denote by rii the cardinality of Jj. 

One will assume that K C M" in ([2]) is compact, and one knows some M > such 
that X G K ^ M - ||xf > 0. For every i < N, introduce the quadratic polynomial 
X I— gm,+i{^) = M — X^fcg/i-^fc- '^^^ index set {1, . . . ,m + N} has a partition U^^Jj 
with Ji ^ ^ for every i = 1, . . . , N. In the sequel we assume that for every i = 1, . . . ,N, 
Pi, Qi G M[xfe : k E li] and for every j G Jj, gj G M[xa; : k E li]. Next, for every i = 1, . . . , N, 
let 

Ki := {z G M"' : (^^(z) > 0, A; G J,} 
so that K in ([2]) has the equivalent characterization 

K = {x G M'' : {xk, k e li) EKi, i = l,..., N}. 

Similarly, for every i,j G {1, . . . , N} such that i j and fl Ij ^ 0, 

= Kji := {(xfc, k e liH Ij) : (xfc, k E U) E K^; (x^. A; G /j) G }. 

Let A^(K) be the space of finite Borel measures on K, and for every i = 1, . . . , A^, let 
TTj : A^(K) — )■ A^(Kj) denote the projection on Kj, that is, for every /i G A^(K): 

7r,/i(5) := /i({x : X G K; (x^. A; G J^) G 5}), V5 G i3(K,) 

where -B(Kj) is the usual Borel cx-algebra associated with Kj. 

For every i,j E {1, . . . , A^} such that i ^ j and Jj fl Ij ^ 0, the projection VTjj : A^(Ki) — )■ 
A^(Kjj) is also defined in an obvious similar manner. For every i = 1, . . . , — 1 define 
the set: 

:= {3 E{i + l,...,N} : /,n/,^0}, 
and consider the infinite dimensional problem 

N 



f := inf / pidfii 



s.t. / Qidni = 1, i = l,...,N i"^) 



7iij{qidfj.i) = TCji{qjdfij), Wj E Ui, i = 1, . . . , N - 1. 



Definition 3.1 Sparsity pattern satisfies the running intersection property if for 

every i = 2, . . . , N : 



Pi ( U 4 j ^ Ij, for some j <i-l. 

\k=l 
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Theorem 3.1 Let K C M" in Ij^ be compact. If the sparsity pattern satisfies the 

running intersection property then f = f* ■ 

Proof: That / < /* is straightforward. As K is compact and > on K for every 
i = 1,...,N, f* = J2iLi fii'^*) some x G K. So let fi be the Dirac measure 6y_* 
at X* and let z/j be the projection TTj/i of /x on Kj. That is z/j = 5i^xi,k&ii)i the Dirac 
measure at the point (x^, k G U) of Kj. Next, for every i = 1, . . . , iV, define the measure 
:= gi(x*)~^ciz/i. Obviously, (/ij) is a feasible solution of ([7j) because iii G A4(Kj) and 
J giC?/Xj = 1, for every i = 1, . . . , A^, and one also has: 

(x^, k E liH Ij) = VTjj/ij = Hji^j, Vj 7^ i such that Ij H li 0. 

Finally, its value satisfies 

N „ N 



/ pidf^i = J^pi(x*)/gi(x*) = /*, 

i=i i=i 



and so / < /*. 

We next prove the converse inequality f > f*. Let (/ij) be an arbitrary feasible solution of 
([7j), and for every i = 1, . . . , N, denote by z/j the probability measure on Kj with density 
qi with respect to /ij, that is. 



(B) := / q,{^)dfi,{^), V5Gi3(K, 
JKinB 



By definition of the linear program ([7]), TTjjZ/j = TTjiiyj for every couple j ^ i such that 
/j n /j 7^ 0. Therefore, by [9j, Lemma B.13] there exists a probability measure z^ on K 
such that TTjZ^ = z/j for every i = 1, . . . , N. But then 



N ^ N ^ N 



/ Pidfii = y^ —dui = / — ciz/ 
and so / > /*. □ 



3.2 A hierarchy of sparse SDP relaxations 

Let y = (ya) be a real sequence indexed in the canonical basis (x") of ]R[x]. Define the 
linear functional Ly : R[x] — )■ M, by: 

For every i = 1, . . . , A^, let 

N(^) := { a G : = if ^ Jj }; nJ!^ := { a G M^^^ : J] ttj < A; }. 
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An obvious similar definition of N'-*-'-' (= N^-^*-*) and n[*'''' (= N^^*"*) applies when considering 

/, n /, ^ 0. 

Let y = (ya) be a given sequence indexed in the canonical basis of M[x]. For every 
i = 1,. . . ,N, the sparse moment matrix Mk{y,Ii) associated with y, has its rows and 
columns indexed in the canonical basis (x"") of M.[xk '■ A; G /«], and with entries: 

Mfc(y, J,)„,^ = i:y(x-+'^) = y^+p, \/a,f3e 

Similarly, for a given polynomial h G : G /j], the sparse localizing matrix 

Mfc(/iy, Jj) associated with y and h, has its rows and columns indexed in the canoni- 
cal basis (x") of M.[xk '■ A; G /j], and with entries: 

With K C M" defined in ([2]), let rj := [(deg(7j)/2] , for every j = 1, . . . ,m + N. Consider 
the hierarchy of semidefinite relaxations: 

N 

s.t. Mfe(y,/i)^0, 

Mfc_,^.((7,y,/,) ^0, 

Ly(x"g,) = Ly(x"g,) 



^ = l,...,iV 

VjG J„2 = l,...,iV (8) 
1 = 1,. ..,N 
= 0, Va G N(^J'), Vj G f/i, i = 1, . . . , - 1 
with |q;| + max[deggj, degg,] < 2k. 



Theorem 3.2 Le^ K C m ([^ 6e compact. Let the sparsity pattern {Ii)fLi satisfy 
the running intersection property, and consider the hierarchy of semidefinite relaxations 
defined in Then: 

(a) fk t f* ask ^oo. 

(b) // an optimal solution y* of ^ satisfies 

rankMfc(y*, /,) = rankMfc_,^(y*, J^) =: R„ Wi = l,...,N, 
(where Vi = maxjg jjr^]^, and 

rankMfc(y*,/,n/,) = 1, Vj G f/„ =l,...,iV-l, 
then fl = f* and one may extract finitely many global minimizers. 

Proof: The proof is similar to that of Theorem 12.21 and also to that of j9| Theorem 4.7]. 
One first prove that if (y'^) is a nearly optimal solution of ([5]), i.e. 

^ 1 
fl<Y.L^^{Pi)<fl + ^i 

1=1 

then there exists a subsequence (kf) and a sequence y, such that 

lim 1/^^ = y,, VaGN«,^ = l,...,A^. 
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From this pointwise convergence it easily follows that for every i = 1, . . . , N and j E Ji. 
Mk{y, h) h 0, Mkigj y, /,) ^0, j e Jf, t = 1, . . . , N. 



Now observe that each set K, C M"' satisfies Assumption 12.11 Therefore, by Putinar's 
theorem [9l Theorem 2.14] the sequence y* = (ya), ce G N^*^ (a subsequence of y), has a 
representing measure /ij supported on Kj. For every with j G Ui, denote by y*-' the 
sequence (ya), ol G N^*-'^ Again, by pointwise convergence, Ly{qi) = 1, i = 1, . . . , N, and 



Ly(gix") = [ x"gj(x)d/ij- = [ x"gj(x)d/ij-, Va G Vj G Ui. 



(9) 



Therefore, for every i = 1,. . . ,N, dvi := qi{'x)dfii is a finite Borel probability measure 
supported on Kj. As measures on compact sets are moment determinate, ([9]) yields: 

TlijUi = TljiUj, V(2, j), j G Ui. 

Therefore, by j9l Lemma B.13] there exists a probability measure z/ on K such that 
TiiU = Ui for every i = 1, . . . , N. But then 

r>lim/4 = lim Vl = V / p^rf/ii 

i=l i=l "'^ 

i=i 1' i=i 



As the converging subsequence was arbitrary, and (/^) is monotone non decreasing, we 
finally get t /*• addition, z/ is an optimal solution of ([7]) with optimal value f* = f- 

The proof of (b) is as in |7] and uses the fiat extension theorem of Curto and Fialkow P 
Theorem 3.7] from which, the sequence y* = (?/*), a G N*-*\ has an atomic representing 
measure supported on Ri points of K,, for every i = 1, . . . , N . □ 



4 GloptiPoly and examples 

In this section we show that the generalized moment problem (GMP) formulation of 
rational optimization problem ([1]) has a straightforward Matlab implementation when 
using our software GloptiPoly 3 [S] . Rather than explaining the approach in full generality 
with awkward notations, we describe three simple examples. 



4.1 Wilkinson-like rational function 

Consider the elementary univariate rational optimization problem 

r = sup/(x), m = f2^ = f2^ 
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with an integer. The only real critical point is x = 0, at which the objective function 
takes its maximum 

^ 1 

r = m = T.f 

i=l 

Reducing to the same denominator 

yields the well-known Wilkinson polynomial q whose squared root moduli are the integers 
from 1 to N. This polynomial was described in the mid 1960s by James H. Wilkinson 
to illustrate the difficulty of finding numerically the roots of polynomials. If we choose 
e.g. N = 20, reduction to the same denominator is hopeless since the constant coefficient 
in monic polynomial q is 20! = 2432902008176640000. The GMP formulation (g]) of this 
problem reads (up to replacing inf with sup in the objective function): 

s.t. /j^gid/ii = 1 

J^x'^qidni = J^x^qidfii, Va G N", i = 1,. . .,N. 

Our Matlab script to model and solve this problem is as follows: 

N = 20; mpol('x',N); % create variables 
q = cell(N,l); % problem data 
mu = cell(N,l); /(, measures 

for i = 1:N, q{i} = i+x(i)"2; mu{i} = meas(x(i)); end 
y„ model GMP 

k = 0; % relaxation order 
f = mass(mu{l}); % objective function 
e = [mom(q{l}) == 1] ; % moment contraints 
for i = 2:N 
f = f + mass(mu{i}); 

e = [e; mom(mmon(x(l) ,k) *q{l}) == mom(mmon(x(i) ,k) *q-[i}-)] ; 
end 

% model SDP relaxation of GMP 
P = msdp(max(f ) , e) ; 
"/o solve SDP relaxation 
[statjObj] = msol(P) 

Instructions mpol, meas, mass, mom, mmon, msdp, max and msol are GloptiPoly 3 commands, 
see the user's guide [5j for more information. For readers who are not familiar with this 
package, variable f is the objective function to be maximized. Since Pi = 1 for all 
z = 1, . . . , A^, it is the sum of masses of measures /ij. Vector e stores the linear moment 
constraints and the instruction mmon(x,k) generates all monomials of variable x up to 
degree k. Finally, instruction msdp generates the SDP relaxation of the GMP, and msol 
solves the SDP problem with the default conic solver (SeDuMi 1.3 in our case). 

At the first SDP relaxation (i.e. for k=0) we obtain a rank-one moment matrix corre- 
sponding to a Dirac at = 0: 

» [stat,obj] = msol(P) 
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Global optimality certified numerically 
Stat = 
1 

obj = 

3.5977 

which is consistent with Maple's 

> f := sum(l/(x"2+i) , i=1..20); 

> evalf(subs(x = 0, f)); 

3.5977 

Note that for this example Assumption 12.11 is violated, since we optimize over the non- 
bounded set K = M. In spite of this, we could solve the problem globally. 



4.2 Relevance of the compactness assumption 

With this elementary example we would like to emphasize the practical relevance of As- 
sumption UTT] on the existence of an algebraic certificate of compactness of set K. Consider 
the univariate problem 

r=inf/(x), f{x)= l + + 2±^. (10) 

First let K = M. The numerator of the gradient of f{x) has two real roots, one of which 
being the global minimum located at x* = —1.4215 for which /* = 1.1286. The following 
GloptiPoly script models and solves the SDP relaxations of orders k = 0, . . . , 9 of the 
GMP formulation of this problem: 



mpol xl x2 

fl = l+xl+xl"2; gl = l+xl"2; f2 = l+x2-2; g2 = l+2*x2-2; 
mul = meas(xl); mu2 = meas(x2); 
bounds = [] ; 
for k = 0:9 

P = msdp(min(mom(f l)+mom(f 2)) , ... 

moni(nimon(xl ,k) *gl) == moni(nimon(x2,k)*g2) , mom(gl) == 1); 

[stat, obj] = msol(P); 

bounds = [bounds ; obj ] ; 
end 

bounds 



In vector bounds we retrieve the following monotically increasing sequence of lower bounds 
fk (^P to 5 digits) obtained by solving the SDP relaxations (jS]): At SDP relaxation k = 9, 
GloptiPoly certifies global optimality and extracts the global minimizer. Table [T] shows 
that the convergence of the hierarchy of SDP relaxations is rather slow for this very simple 
example. This is due to the fact that Assumption 12.11 is violated, since we optimize over 
the non-bounded set K = M. 

On Figure [1] we report the sequences of lower bounds obtained by solving the SDP relax- 
ations of problem ffTOj) on compact sets K = [—R, R] for i? = 2, 3, . . . , 9. 
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order k 


bound 


order k 


bound 





1.0000 


5 


1.0793 


1 


1.0000 


6 


1.1264 


2 


1.0170 


7 


1.1283 


3 


1.0220 


8 


1.1286 


4 


1.0633 


9 


1.1286 



Table 1: Lower bounds for SDP relaxations of problem f llOp . 



1.15 



1.1 



1.05 




1 23456789 10 

SDP relaxation order 

Figure 1: Lower bounds for SDP relaxations of problem (ITUl) on bounded sets K = [—R, R] 
for i? = 2 (top curve) to i? = 9 (bottom curve). 

4.3 Exploiting sparsity with GloptiPoly 

Even though version 3 of GloptiPoly is designed to exploit problem sparsity, there is no 
illustration of this feature in the software user's guide [S]. In this section we provide such 
a simple example. Note also that GloptiPoly is not able to detect sparsity in a given 
problem, contrary to SparsePOP which uses a heuristic to find chordal extensions of 
graphs [12]. However, SparsePOP is not designed to handle directly rational optimization 
problems. 

Consider the elementary example of [3 Section 3.2]: 

infxgR4 X1X2 + X1X3 + XiXa 
s.t. x\ + xl<l 
xl + xl<2 
~\~ 3 

for which the variable index subsets Ji = {1,2}, I2 = {1,3}, I3 = {1,4} satisfy the 
running intersection property of Definition 13.11 Note that this problem is a particular 
case of ([1]) with a polynomial objective function. 

Without exploiting sparsity, the GloptiPoly script to solve this problem is as follows: 
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mpol xl x2 x3 x4 

Pdense = msdp(min(xl*x2+xl*x3+xl*x4) , ... 

xl~2+x2"2<=l,xl"2+x3~2<=2,xl"2+x4"2<=3,2) ; 
[statjObj] = msol (Pdense) ; 

GloptiPoly certifies global optimality with a moment matrix of size 15, and 3 localizing 
matrices of size 5. And here is the script exploiting sparsity, splitting the variables into 
several measures fii consistently with subsets If. 

mpol xl 3 
mpol x2 x3 x4 

mu(l) = meas([xl(l) x2] ) ; 7„ first measure on xl and x2 

mu(2) = meas([xl(2) x3] ) ; % second measure on xl and x3 

mu(3) = meas([xl(3) x4] ) ; 7„ third measure on xl and x4 

f = mom(xl(l)*x2)+mom(xl(2)*x3)+mom(xl(3)*x4) ; % objective function 

k = 3; % SDP relaxation order 

ml = mom(mmon(xl (1) ,k) ) ; % moments of first measure 

m2 = mom(mmon(xl (2) ,k) ) ; % moments of second measure 

m3 = mom(mmon(xl (3) ,k) ) ; % moments of third measure 

K = [xl(l)"2+x2~2<=l, xl(2)"2+x3"2<=2, xl (3) "2+x4"2<=3] ; I supports 

Psparse = msdp(min(f ) ,ml==m2,m3==m2,K,mass(mu)==l) ; 

[stat,obj] = msol (Psparse) ; 

GloptiPoly certifies global optimality with 3 moment matrices of size 6, and 3 localizing 
matrices of size 3. 



4.4 Comparison with the epigraph approach 

In most of the examples we have processed, the epigraph approach described in the In- 
troduction (consisting of introducing one lifting variable for each rational term in the 
objective function) was less efficient than the GMP approach. Typically, the order of the 
SDP relaxation (and hence its size) required to certify global optimality is typically larger 
with the epigraph approach. 

When evaluating the epigraph approach, we also observed that it is numerically preferable 
to replace the inequality constraints rjgj(x) — pj(x) > with equality constraints rjgj(x) — 
Pj(x) = in the definition of semi-algebraic set K in ([3]). For the example of Section |4?T] 
the epigraph approach with inequalities certifies global optimality at order k = 5, whereas 
the epigraph approach with equalities requires k = 1. 

As a typical illustration of the issues faced with the epigraph approach consider the 
example with eighth-degree terms 
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mf f(x) = > ^ ^^-^ \^ „ o n . (11) 

^^^^ fr^ {x\ + xl + 2i){xl + xlxl + x^ + i^) ^ ' 

which is cooked up to have several local optima and sufficiently high degree to prevent 
reduction to the same denominator. After a suitable scaling to make critical points fit 
within the box [—1,1]^, as required by the moment SDP relaxations formulated in the 
power basis [H Section 6.5], the GMP approach yields a certificate of global optimality 
with x\ = —0.60450, X2 = —2.2045, /* = —6.2844 at order = 6 in a few seconds on 
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a standard PC. In contrast, the epigraph approach does not provide a certificate for an 
order as high as k = 10, requiring more than one minute of CPU time. 



4.5 Shekel's foxholes 

Consider the modified Shekel foxholes rational function minimization problem [2] 

^ 1 

min > — — (12) 

whose data aij, Ci, i = 1,...,N, j = l,...,n can be found in [U Table 16]. This 
function is designed to have many local minima, and we know the global minimum in 
the case n = 5, see [H Table 17]. After a suitable scaling to make critical points fit 
within the box [0, 1]^, and after addition of a Euclidean ball constraint centered in the 
box, the GMP approach yields a certificate of global optimality at order /c = 3 in less 
than one minute on a standard PC. The extracted minimizer is a;* = 8.0254, = 9.1483, 
Xg = 5.1138, X4 = 7.6213, = 4.5638, which matches with the known global minimizer 
to four significant digits. This point can be refined if given an initial guess for a local 
optimization method. If we use a standard quasi-Newton BFGS algorithm, we obtain 
after a few iterations a point matching the known global minimizer to eight significant 
digits. 

In the case n = 10, for which the global minimum is given in |T1 Table 17], the GMP 
approach yields a certificate of global optimality at order A; = 2 in about 750 seconds 
of CPU time. Here too, we observe that the extracted minimizer = 8.0249, = 
9.1518, X* = 5.1140, X* = 7.6209, x* = 4.5640, x* = 4.7110, x* = 2.996, x* = 6.1259, 
Xg = 0.73424, x*0 = 4.9820 is a good approximation to the minimizer, with four correct 
significant digits. If necessary, this point can be used as an initial guess for refining with 
a local solver. 

Note that it is not possible to exploit problem sparsity in this case, since all the variables 
appear in each term in sum f|T2l) . 



4.6 Rosenbrock's function 

Consider the rational optimization problem 

n— 1 ^ 

f* — j]2ax N (13) 

xgM" ^ 100(Xi+i - X2)2 + (Xi - 1)2 + 1 ^ ^ 

which has the same critical points as the well-known Rosenbrock problem 

n-l 

min J](100(x,+i-x2)2 + (a;,-l)2) 

i=l 

whose geometry is troublesome for local optimization solvers. It can been easily shown 
that the global maximum /* = 1 of problem fll3p is achieved at x* = 1, i = 1, . . . , n. Our 
experiments with local optimization algorithms reveal that standard quasi-Newton solvers 
or functions of the Optimization toolbox for Matlab, called repeatedly with random initial 
guesses, typically yield local maxima quite far from the global maximum. 
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With our GMP approach, after exploiting sparsity and adding bound constraints xf < 16, 
i = 1, . . . ,n, we could solve problem f|T3l) with a certificate of global optimality for n up to 
1000. Typical CPU times range from 10 seconds for n = 100 to 500 seconds for n = 1000. 



5 Conclusion 

The problem of minimizing the sum of many low-degree (typically non-convex) rational 
fractions on a (typically non-convex) semi-algebraic set arises in several important ap- 
plications, and notably in computer vision (triangulation, estimation of the fundamental 
matrix in epipolar geometry) and in systems control {H2 optimal control with a fixed- 
order controller of a linear system subject to parametric uncertainty). These engineering 
problems motivated our work, but the application of our techniques to computer vision 
and systems control will be described elsewhere. These fractional programming problems 
being non convex, local optimization approaches yield only upper bounds on the optimum. 

In this paper we were interested in computing the global minimum (and possibly global 
minimizers) or at least, computing valid lower bounds on the global minimum, for frac- 
tional programs involving a sum with many terms. We have used a semidefinite program- 
ming (SDP) relaxation approach by formulating the rational optimization problem as an 
instance of the generalized moment problem (GMP). In addition, problem structure can 
be sometimes exploited in the case where the number of variables is large but sparsity is 
present. Numerical experiments with our public-domain software GloptiPoly interfaced 
with off-the-shelf semidefinite programming solvers indicate that the approach can solve 
problems that can be challenging for state-of-the-art global optimization algorithms. This 
is consistent with the experiments made in |3j where the (dense) SDP relaxation approach 
was first applied to (polynomial) optimization problems of computer vision. 

For larger and/or ill-conditioned problems, it can happen that GloptiPoly extracts from 
the moment matrix a minimizer which is not very accurate. It can also happen that Glop- 
tiPoly is not able to extract a minimizer, in which case first-order moments approximate 
the minimizer (provided it is unique, which is generically true for rational optimization). 
The approximate minimizer can be then input to any local optimization algorithm as an 
initial guess. 

A comparison of our approach with other techniques of global optimization (reported e.g. 
on Hans Mittelmann's or Arnold Neumaier's webpages) is out of the scope of this paper. 
We believe however that such a comparison would be fair only if no expert tuning is re- 
quired for alternative algorithms. Indeed, when using GloptiPoly the only assumption we 
make is that we know a ball containing the global optimizer. Besides this, our results are 
fully reproducible (Matlab files reproducing our examples are available upon request) and 
our SDP relaxations are solved with general-purpose semidefinite programming solvers. 
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